Análisis de las observaciones asimiladas

files <- Sys.glob("E1/diagfiles/asim*")

obs <- lapply(files, function(f) {
  date_time <- ymd_hms(substr(basename(f), 11, 24))
  ens <- substr(basename(f), 29, 31)
  fread(f, na.strings = c("0.100E+11", "-0.100E+06", "-99999.90", "-100000.00")) %>% 
    .[, date.time := date_time] %>% 
    .[, ens := ens] %>% 
    .[]
}) %>% 
  rbindlist() %>% 
  .[, c("V2", "V4") := NULL]

colnames(obs) <- c("var", "stationID", "type", "dhr", "lat", "lon", "pressure", "usage.flag", "obs", "obs.guess", "obs2", "obs.guess2", "rerr", "date.time", "ens")

Si analizamos la cantidad de observaciones para cada tiempo, no se ve mucha diferencia. Los tipos 131, 181 y 281 son los mas abundantes (seguramente la mayoría de obs 181 y 281 provienen de las estaciones automáticas que agregamos). Pero en realidad acá tenemos los 10 miembros del emsable juntos, si contamos la cantidad de observaciones en el diag de cada miembro nos da lo mismo. Suena raro pero revisé los diag y si bien tienen la misma cantidad de observaciones, la diferencia entre la observación y el guess cambia. Es posible que haya algo en el código que primero seleccione las observaciones a comparar con el guess? Tal vez ese trabajo lo hace cuando calcula el diag para la media.

obs[, .N, by = .(var, type)] %>% 
  ggplot(aes(factor(type), N)) +
  geom_col(position = "dodge") +
  # scale_fill_viridis_d() +
  scale_y_log10() +
  facet_wrap(~var, scales = "free") +

obs[var == "t" & type == 181 & ens == "001", .N, by = .(stationID, lon, lat)] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_sf(data = map, inherit.aes = FALSE) +
  geom_point(aes(fill = N, size = N), shape = 21) +
  scale_size_area(max_size = 3, guide = "none") +
  coord_sf(ylim = c(-42, -20), xlim = c(-75, -52)) +
  labs(title = "Distribución de estaciones y \ncantididad de observaciones en cada una", x = "lon") +
  theme(legend.position = "bottom")

En los diag se guarda la diferencia entre la observación y el guess. Si todo anda más o menos bien sería deseable que los valores de esta diferencia esten centrados alrededor de 0 con una distribucón más o menos gausiana.

obs[usage.flag == 1] %>% 
  ggplot(aes(obs.guess)) +
  geom_density(aes(color = var)) +
  geom_rug() +
  facet_wrap(~factor(var), scales = "free") +
  theme(legend.position = "none") +
obs[usage.flag == 1] %>% 
  ggplot(aes(obs.guess, rerr)) +
  geom_point(alpha = 0.05, size = 0.5) +
  facet_wrap(~var, scales = "free")

Es importante filtrar las observaciones con usage.flag == 1 porque son las ¿asimiladas? (no se si enkf las usa directamente o hace un nuevo control de calidad). Las obs con usage.flag == -1 es claramente es horrrible y sería interesante ver el origen de esas observaciones. Pero mientras, la distribución parece más o menos razonable aunque con bias positivo o negativo dependiendo de la variable. Es preocupante la asimeria de la temperatura, donde se ve una diferencia entre la observación y el guess de 20 o 30 grados!

Los diags también guardan una medida del error, de donde sale eso? Es un poco más difícil de responder. En principio toma el error de la observación que vienen en el prepbufr (caso GLOBAL) o de la tabla errtable (caso REGIONAL). Pero luego define un nuevo error buscando el máximo entre el error original y el error de la tabla convinfo ¿que va a ser utilizando en el gross check?. Lo que guarda en el diag es el error original o si es muy si es muy pequeño un huge_single que es básicamente cero.

No tengo idea que está pasando con las observaciones de temperatura que tienen error 0 pero son los que tienen una diferencia mayor con el guess. Vamos a mirarlas más de cerca analizando de donde vienen.

181 es el tipo de observaciones de temperatura y humedad que agrego manualmente al prepbufr. Es posible que no todas las observaciones 181 sean de la red de estaciones automáticas porque no revisé el prepbufr antes de agregarlas, pero seguro son la mayoría. Y estas observaciones las causantes de la asimetría en la distribución del obs.guess de la temperatura; aunque definitivamente también se ven cosas raras para las observaciones 180!

Al revisar la tabla errtable se ve que las observaciones 181 tienen un error de 1.5ºC para todos los niveles, con lo cual las observaciones con error cero que se ve, seguramente no cumplen alguno de los requisitos y le asignan un nuevo error. De todas maneras no parece ser un problema constante.

obs[usage.flag == 1 & var == "t" & abs(dhr) <= 0.5] %>% 
  ggplot(aes(obs.guess)) +
  geom_density(aes(group = type, color = factor(type)))

Vamos a aislar las observaciones de temperatura que caen en la cola de la distribución.

Las casi 700 observaciones “malas” están distribuidas en solo 13 estaciones (de más de 400). Hay una “INT511” que tiene más de 150 de estas observaciones malas!

stations <- obs[usage.flag == 1 & var == "t" & type == 181 & obs.guess > 10] %>% 
  group_by(stationID) %>% 
  filter(row_number(stationID) == 1)

obs[usage.flag == 1 & var == "t" & type == 181 & obs.guess > 10 & abs(dhr) <= 0.5, .N, by = .(stationID, lon, lat)] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_sf(data = map, inherit.aes = FALSE) +
  geom_point(aes(size = N), alpha = 0.5) +
  ggrepel::geom_label_repel(data = stations, aes(label = stationID), nudge_y = 2, alpha = 0.5) +
  labs(x = "lon") +
  coord_sf(ylim = c(-42, -20), xlim = c(-75, -52)) 

obs[usage.flag == 1, .(rmsi = sqrt(sum(obs.guess)^2/.N)), by = .(lon, lat, date.time)] %>% 
  copy() %>% 
  .[, date.time := factor(date.time)] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_sf(data = map, inherit.aes = FALSE) +
  geom_point(aes(color = rmsi)) +
  scale_color_distiller(name = "RMSI", type = "seq", 
                        palette = "YlGnBu", direction = -1,
                        limits = c(0, 50)) +
  coord_sf(ylim = c(-42, -20), xlim = c(-75, -52)) +
  facet_wrap(~ date.time) +
  theme(legend.position = "bottom")

Analisis de los análisis y el emsamble

files <- Sys.glob("E1/ANA/*/analysis.mem*")

ana_T2 <- lapply(files, function(f) {
  date_time <- ymd_hms(substr(f, 8, 21))
  ens <- substr(basename(f), 13, 15)
  ReadNetCDF(f, vars = c("T2", "XLAT", "XLONG")) %>% 
    .[, c("Time", "south_north", "west_east") := NULL] %>% 
    .[, date.time := date_time] %>% 
    .[, ens := ens] %>% 
    .[]
}) %>% 
  rbindlist()

colnames(ana_T2) <- c("var", "lat", "lon", "date.time", "ens")

Calculemos la media y el spread del ensamble para cada tiempo.

ana_T2_sum <- ana_T2[, .(mean = mean(var), spread = sd(var)), by = .(lon, lat, date.time)]

ana_T2_sum[mean > 260] %>% 
  copy() %>% 
  .[, date.time := factor(date.time)] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = mean - 273.15)) +
  scale_color_distiller(name = "Mean", type = "seq", palette = "YlGnBu", direction = -1) +
  geom_sf(data = map, inherit.aes = FALSE, fill = NA) +
  coord_sf(ylim = c(-42, -20), xlim = c(-75, -52)) +
  facet_wrap(~ date.time) +
  theme(legend.position = "bottom") +

ana_T2_sum[mean > 260 & spread < 2.5] %>% 
  copy() %>% 
  .[, date.time := factor(date.time)] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = spread)) +
  scale_color_distiller(name = "Spread", type = "seq", palette = "YlGnBu", direction = -1) +
  geom_sf(data = map, inherit.aes = FALSE, fill = NA) +
  coord_sf(ylim = c(-42, -20), xlim = c(-75, -52)) +
  facet_wrap(~ date.time) +
  theme(legend.position = "bottom")

Hay algunos poquitos puntos que se van de escala con temperaturas de casi -30º y spread super alto. Los saco para ver el resto de los detalles pero habría que revisar que puede estar pasando. Si analizamos la diferencia entre la media del ensamble y cada miembro –> anom también encuentro que hay algunos puntos con valores muy altos y eso no me deja ver mucho.

ana_T2[, c("mean", "spread") := .(mean(var), sd(var)), by = .(lon, lat, date.time)] %>% 
  .[, anom := var - mean] %>% 
  copy() %>% 
  .[, date.time := factor(date.time)] %>% 
  .[abs(anom) < 3] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = anom)) +
  scale_color_divergent(name = "Anomalia") +
  geom_sf(data = map, inherit.aes = FALSE, fill = NA) +
  coord_sf(ylim = c(-42, -20), xlim = c(-75, -52)) +
  facet_grid(date.time ~ ens) +
  theme(legend.position = "bottom")

Alguien me explica por favor porque los valores altos de anomalía están en la costa???!!!

ana_T2[, c("mean", "spread") := .(mean(var), sd(var)), by = .(lon, lat, date.time)] %>% 
  .[, anom := var - mean] %>% 
  copy() %>% 
  .[, date.time := factor(date.time)] %>% 
  .[abs(anom) > 3] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = anom)) +
  scale_color_divergent(name = "Anomalia") +
  geom_sf(data = map, inherit.aes = FALSE, fill = NA) +
  coord_sf(ylim = c(-42, -20), xlim = c(-75, -52)) +
  facet_wrap(~ens, ncol = 5) +
  theme(legend.position = "bottom")